# Function to prepare Data and run stan to sample from posteror 
######################

est_ros_stan <- function(data,iteration=200, nchains=1,nburnin=100, country, year,
                     pos="lr",
                     model="model_rozenas_graded_response.stan",
                     ...){
  
  
  # Install required Packages
  require(reshape2)
  require(rstan)  
  rstan_options(auto_write = TRUE) # Multipel Cores
  options(mc.cores = nchains)
  
  # Prepare Data
    sel  <- c("ctry_name","pid","eid",pos)
    data <- data[data$ctry_name==country & data$year==year,]
    
    # Drop Experts that place all as missing
    ld_eid <- split(data, data$eid)
    am_eid <- sapply(ld_eid, function(d)  all(is.na(d$lr)) )
    sel_row <- !(data$eid %in% names(am_eid)[am_eid])
    data <- data[sel_row,]  
    
  # Position Data
    cd   <- melt(data[,sel], id.vars = c("ctry_name","pid","eid"))
    cd   <- na.omit(cd)
    
    # recode items 1:K
    cd$variable <- as.character(cd$variable)
    ui <- unique(cd$variable)
    for(i in 1:length(ui)) cd$variable[cd$variable==ui[i]] <- i
  
  # Missigness Data
    data$pos_mis <- as.numeric(is.na(data[,pos]))
    sel  <- c("ctry_name","pid","eid","pos_mis")
    cdm <- melt(data[,sel], id.vars = c("ctry_name","pid","eid"))
  
    
  # Stan Data
    uc <- function(x) length(unique(as.factor(x)))
    standata <- list(T = 11, 
                     J = uc(cd$pid),
                     N = nrow(cd), 
                     N_miss = nrow(cdm),
                     I = uc(cd$eid),
                     jj = as.numeric(as.factor(cd$pid)),
                     ii = as.numeric(as.factor(cd$eid)),
                     jj_miss = as.numeric(as.factor(cdm$pid)),
                     ii_miss = as.numeric(as.factor(cdm$eid)),
                     y = cd$value +1,
                     y_miss = cdm$value,
                     c = seq(from=-1,to=1,length=11-1)
                     )
  
  
  # Initial Value Function; lstoetze: nu is fixed otherwise nummerical errors
  init_fun <- function(T = uc(cd$value), J = uc(cd$pid), K = uc(cd$variable),I = uc(cd$eid)) {
    list(
       mu = rnorm(J, 0.1)
      ,nu = rep(1,J)
      ,phi  = rep(1,I)
      ,tau = rnorm(I,0,.1)
      ,alpha = rnorm(2,0,0.1)
    )
  }

  # Call stan
  fit <- stan(file = model, data = standata
              ,init = init_fun
              ,iter = iteration, chains = nchains,warmup = nburnin,
              ...)
  
   # Return Fit
  return(fit)
}

